Mathematical modeling of circadian rhythms#
circStudio implements four light-informed mathematical models designed to characterize and predict circadian physiology: Forger, Jewett, HannaySP and HannayTP. All four models can be used to estimate core body temperature minimum and the dim light melatonin onset. We begin by importing the necessary libraries:
import circStudio
import os
Import sample RPX file:
# Create file path for sample files within the circStudio package
fpath = os.path.join(os.path.dirname(circStudio.__file__))
# Create a new Raw instance using awd RPX adaptor
raw = circStudio.io.read_rpx(os.path.join(fpath, 'data', 'test_sample_rpx_eng.csv'),
start_time='2015-07-04 12:00:00',
light_mode='White Light',
period='6 days',
language='ENG_UK')
1. Dynamic modeling of actigraphy-based light intensity data#
After loading the sample data, the next step is to create an instance of a circadian model for analysis. The model can be initialized using either a pd.Series or light data indexed by datetime, or two separate np.array objects, one for time and one for light intensity.
# Create a new instance of the HannaySP model using light data
model = circStudio.HannaySP(data=raw.light)
In this tutorial, we will demonstrate the functionality mainly using the HannaySP model, but syntax is consistent across the other available models:
# Jewett model
model = circStudio.Jewett(data=raw.light)
# Forger model
model = circStudio.Forger(data=raw.light)
# HannayTP model
model = circStudio.HannayTP(data=raw.light)
The resulting model states are stored in the model_states attribute. Be default, when a model is instantiated, the trajectory is computed using a predefined set of initial conditions, assuming a light schedule of 16 hours of light followed by 8 hours of complete darkness. To obtain more realistic and representative initial conditions for the model states, the model can be run iteratively over the light vector util it reaches a stable solution, assessed by convergence of the DLMO. For this, users must specify the number of iterations (loop_number) and provide the light data to be used in the loop:
# Loop over the provided light trajectory to determine new initial state
model.get_initial_conditions(data=raw.light, loop_number=50)
The model entrained after 6 loops.
array([ 0.77016547, 263.37108225, 0.26764411])
From the model output, it is possible to extract a series of predicted DLMOs. By default, these values are expressed in hours relative to the start of the recording. To convert them to a 24-hour clock format, the modulo operator (% 24) can be applied:
# Obtain array of daily DLMOs
model.dlmos() % 24
array([7.525 , 7.4 , 7.28333333, 6.98333333, 6.7 ,
7.3 ])
Using the plot() method in the Model class, it is possible to visually inspect the evolution of the daily predicted DLMO:
model.plot(dlmo=True)
Using the same function, it is possible to visualize all model states:
model.plot(states=True)
A key advantage of the HannaySP model, compared to the other available models, is its ability to compute the Entrainment Signal Regularity Index (ESRI). This metric quantifies how effectively a given light input can entrain the circadian system, ranging from 0 (no entrainment) to 1 (a maximally entraining signal).
# Compute ESRI
esri = circStudio.ESRI(data=raw.light)
# ESRI mean
esri.mean
np.float64(0.595468534541619)
# ESRI std
esri.std
np.float64(0.054790277969275325)
Future versions of circStudio will include additional mathematical models of circadian rhythms, as well as a model comparison utility. This functionality will allow users to systematically compare Kronauer-based, Van der Pol-type models with the physiologically-informed Hannay formalism. Stay tuned for upcoming developments!
2. Cosinor analysis#
# Create a cosinor object
cosinor = circStudio.Cosinor()
# Fit cosinor to raw activity data
results = cosinor.fit(data=raw.activity, verbose=False)
You can access individual best fit parameter values using the following syntax:
# Mesor
mesor = results.params['Mesor'].value
# Acrophase
acrophase = results.params['Acrophase'].value
# Period
period = results.params['Period'].value
# Amplitude
amplitude = results.params['Amplitude'].value
# Goodness-of-fit metrics (Akaike information criterium and Reduced Chi^2)
aic = results.aic
redchi = results.redchi
Additionally, you can also retrieve best fit parameters as a dictionary:
results.params.valuesdict()
{'Amplitude': np.float64(30.32292038752961),
'Acrophase': np.float64(3.6510311570819654),
'Period': np.float64(1626.5421520119878),
'Mesor': np.float64(168.64347298661386)}
Finally, you can generate an interactive plot using the plot() method from Cosinor class:
cosinor.plot(data=raw.activity, best_params=results.params)
3. Next steps#
In the current version of circStudio (v1.0), support for Functional Linear Modeling (FLM), (Multifractal) Detrended Fluctuation Analysis (MF-DFA), and Singular Spectrum Analysis (SSA) has been temporarily removed to allow for reimplementation and enhanced functionality. These methods will be reintroduced in future versions of circStudio. Again, stay tuned!